We are IntechOpen, 
the world’s leading publisher of 


Open Access books 
Built by scientists, for scientists 


5.300 130,000 150M 


ailable International authors and editors Downloads 


Our author among the 


154 TOP 1% 12.2% 


Countries delivered to most cited s Contributors from top 500 universities 


Selection of our books indexed in the Book Citation Index 
in Web of Science™ Core Collection (BKCI) 


Interested in publishing with us? 
Contact book.department@intechopen.com 


Numbers displayed above are based on latest data collected. 
F information visit www.intechopen.com 


ay 


Chapter 12 


Topological Characterization and Advanced Noise- 
Filtering Techniques for Phase Unwrapping of 
Interferometric Data Stacks 


Pasquale Imperatore and Antonio Pepe 


Additional information is available at the end of the chapter 


http://dx.doi.org/10.5772/61847 


Abstract 


This chapter addresses the problem of phase unwrapping interferometric data stacks, ob- 
tained by multiple SAR acquisitions over the same area on the ground, with a twofold 
objective. First, a rigorous gradient-based formulation for the multichannel phase un- 
wrapping (MCh-PhU) problem is systematically established, thus capturing the intrinsic 
topological character of the problem. The presented mathematical formulation is consis- 
tent with the theoretical foundation of the discrete calculus. Then within the considered 
theoretical framework, we formally describe an innovative procedure for the noise filter- 
ing of time-redundant multichannel multilook interferograms. The strategy underlying 
the adopted multichannel noise filtering (MCh-NF) procedure arises from the key obser- 
vation that multilook interferograms are not fully time consistent due to multilook opera- 
tions independently applied on each single interferogram. Accordingly, the presented 
MCh-NF procedure suitably exploits the temporal mutual relationships of the interfero- 
grams. Finally, we present some experimental results on real data and show the effective- 
ness of our approach applied within the well-known small baseline subset (SBAS) 
processing chain, thus finally retrieving the relevant Earth’s surface deformation time ser- 
ies for geospatial phenomena analysis and understanding. 


Keywords: SAR interferometry, phase unwrapping, discrete calculus 


1. Introduction 


Multichannel (or multitemporal) InSAR techniques address the processing of interferometric 
data stacks obtained by combining multiple SAR acquisitions over the same area. These 
approaches can be essentially categorized in two main classes: persistent scatterers (PS) and 
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small baseline (SB)-based techniques. The solution of the multichannel phase unwrapping 
(MCh-PhU) problem is generally required in multichannel InSAR techniques, whenever 
multidimensional SAR data set, conveying information about complex Earth’s crust events, 
have to be systematically investigated on suitable space-time scales for geospatial phenomena 
understanding [1-29]. In this chapter, we focus on two different related main issues. 


Primarily, we present a rigorous gradient-based formulation of the MCh-PhU problem that is 
consistent with the theoretical foundation of the discrete calculus [30-34]. Emphasis is placed 
on the topological characterization of the underlying discrete setting provided by the differen- 
tial operators of the discrete calculus, which are formally associated with matrix operators and 
represent the discrete counterparts of the classical differential operators of the vector calculus. 
Accordingly, MCh-PhU problem formulation is systematically established in terms of discrete 
differential operators, which are defined by the topology of the intrinsically discrete spaces 
upon which they act, thus capturing the essential topological character of the problem within 
a systematic matrix formalism [35]. It is worth highlighting that our approach provides an 
unambiguous and theoretical-consistent formalism for the MCh-PhU problem, overcoming 
the conceptual inconsistencies of the existing gradient-based formulations [1, 17, 29]. Indeed, 
the existing approaches pose some conceptual limitations from a mathematical viewpoint since 
they rely on an intrinsically discrete setting based description and, at the same time, resort to 
the concepts of gradient and curl of the conventional vectorial calculus, which inherently imply 
a reference to an underlying continuum space and the notion of the infinitesimal [30]. In 
addition, the proposed formal framework enables meaningful analytical investigations on a 
mathematical consistent playground, also providing interesting implications and permitting 
to express previous obtained results in a more general way. 


Then we present an innovative procedure to filter out the noise affecting the phase components 
of a redundant set of (multitemporal) multilook small-baseline interferograms. This is 
achieved by independently solving, for each pixel of the scene, a nonlinear optimization problem 
based on computing the wrapped phase vector that minimizes the (weighted) circular variance 
of the difference between the original and noise-filtered interferograms [43]. This noise- 
filtering procedure arises from the key observation that multilook interferograms are not fully 
time consistent because they are generated through multilook operations that are independ- 
ently carried out on each single interferogram. Indeed, the wrapped discrete curl of the 
interferometric phases defined on a graph whose nodes and edges describe SAR acquisitions 
(in the time/perpendicular baseline domain) and inherent interferograms, respectively, is 
different from zero. This modulo-27 cyclic inconsistency of multichannel interferometric 
phases is properly handled by the presented multichannel noise-filtering (MCh-NF) proce- 
dure. The presented technique is very easy to implement because it does not imply a prelimi- 
nary time-consuming selection of statistically homogenous pixels (SHP), as for instance 
required by the SqueeSAR technique [44], and it has no need of any a priori information on the 
statistics of complex-valued SAR images. The effectiveness of the presented noise-filtering 
approach as well as its impact on the quality of multichannel phase unwrapping procedures 
are also fully investigated. 
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2. Multichannel phase unwrapping problem 


In this Section, we review the mathematical formulation of the multichannel phase unwrap- 
ping (MCh-PhU) problem within the purview of the discrete calculus. As a matter of fact, a 
graph-based description naturally arises in formulating the MCh-PhU problem due to the 
underlying discrete irregular data structure. Indeed, as far as discrete settings (e.g., graphs) 
are concerned, resorting to the conventional vectorial calculus might not be adequate since it 
inherently implies a reference to an underlying continuum space. On the contrary, discrete 
calculus offers a rigorous methodological framework since it treats a discrete domain as entirely 
its own entity. In particular, discrete calculus provides proper differential operators that make 
it possible to purely operate onto a finite, discrete structure without referring to the continuous 
space and notion of the infinitesimal [30]. More specifically, the introduction of some well- 
known algebraic structures [30-33] capturing the essential topological character of the 
underlying graphs permits to phrase pertinent differential operators as matrices. Therefore, one 
of the most important consequences of this approach is that the purely topological nature of 
the discrete differential operators is made more apparent and concrete. Accordingly, by 
systematically adopting the relevant key concepts and tools available within this theoretical 
framework, we here provide a description of the MCh-PhU problem on a suitable mathemat- 
ical playground. For such a purpose, we first establish the notation and terminology used 
throughout the subsequent Sections 2.1, and then the problem at hand is reviewed within this 
formalism (Sections 2.2 and 2.3). We remark that the focus here is on presenting key concepts 
that are useful for the following analyses; however, a comprehensive treatment of the discrete 
calculus and related huge fields of mathematics (e.g., algebraic topology, exterior calculus, and 
differential forms) is clearly beyond the scope of this work but can be found in refs. [30-33]. 


2.1. The theoretical framework of discrete calculus 


A graph GV) is defined by two sets: Y and €. The former is the set of nodes (or vertices) of 
the graph, and the latter represents the corresponding set of edges. Let Q and M be the 
cardinality of V and £, respectively. The vector space R” is referred to as the edge space, and 
the vector space RS is referred to as the vertex space, with R denoting the field of real numbers. 
Without loss of generality, we here assume that the graph @ is connected (i.e., every pair of 
vertices in the graph is connected [33]). Moreover, an orientation establishes a default direction 
on an edge that is considered positive or negative, thus yielding an oriented graph. The M x 
Q incidence matrix II=(II,,,] of an oriented graph g specifies its edge—node connectivity 
relations, whose entries are defined as follows [30-33]: 


-1, if mth edge starts at gth node 
IT, =\+1, if mth edge ends at qth node (1) 
0, otherwise 


with m = 1, 2,..., M and q =1, 2,..., Q. It is important to note that the column rank of IT is Q - 1. 
The incidence matrix IT generates an orthogonal decomposition R(T) e (IT')=R™, where 
R(IT) is the column space of the incidence matrix IT, and MII") denotes the kernel (or null- 
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space) of the matrix I’. The notion of cycle space is also fundamental in graph theory. The 
cycle space of the graph G namely, C= AA, is the subspace of edge space R“ spanned by all the 
cycles in g. The dimension of (( is also referred to as the cyclomatic number of G [33]. It is also 
well known that, for every connected graph G with Q nodes and M edges, the dimension of the 
cycle space is given by R = dim(C()) =M - Q + 1 [30]. Each basis for C(g, (i.e., the cycle basis) 
is therefore uniquely specified by an M x R matrix Q, called cycle matrix. Thus, the column 
vectors of Q=[w!, ..., w*] form a basis for an R-dimensional vector subspace (the cycle space 
of (A of R“. CA is indeed the null-space of I" so that a cycle basis provides a basis for MI) 
[30]. Accordingly, a fundamental property of a linear graph is expressed by the remarkable 
relations: 


'Q=0 (2) 


Q' 11 =0 (3) 


Indeed, several methods for defining a cycle set have been studied, and they can be used to 
define incidence relations between edges and cycles. Specifically, the definition of a cycle set 
from the edge set can be obtained algebraically and geometrically (i.e., from an embedding). 
Algebraic methods find a suitable M x R matrix Q whose columns provide a basis for the null- 
space of IT", with R=dim( MUI" ). Geometric methods for defining a cycle set (i.e., from an 
embedding) permit to identify algorithmically a cycle set (representing the faces) in this 
embedding and may be used to produce the edge—face incidence matrix Q (as illustrated in 
Figure 1). In particular, it is possible to consider a normalized irreducible cycle basis forming 
elementary (or irreducible) cycles [30-33], i.e., cycles that contain no other cycles, so that we 
can associate to each elementary cycle an elementary cycle vector w =[w4 ..., @y]', whose 


entries are defined as follows: 


+1 if rth cycle traverses ith edge forward 


@,=5-1 ifrth cycletraverses ith edge backward (4) 
0 otherwise 
Accordingly, the so defined Q=[w!, ..., w] provides a particularly attractive basis for 


MT"), i.e, the cycle basis formed by all the elementary cycle vectors associated with the 
elementary cycles in g. We also note that Q defines the incidence connectivity relations 
between edges and cycles (see Figure 1). 


It is instructive to highlight that the topological operators I, IT", and Q provide the discrete 
counterparts of the classical gradient (V ), divergence (V -),and curl (V x ) operators of the vector 
calculus for continuous space, respectively. Accordingly, they can be regarded as differential 
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operators on the discrete setting [30]. In addition, it is worth emphasizing that identities (2) 
and (3) mimic the properties of their classical vector calculus counterparts V - V x =0 (div 
curl = 0) and V xV =0 (curl grad = 0), respectively. It should be pointed out that I yields 
differences along edges of nodal “potentials” represented by IIx. Conversely, given an 
arbitrary f € R”, a solution of the equation Ix=f (if it exists) is called the potential of f. Note 
also that x (if it exists) is not unique since the constant column I=(1, ..., 1]'€ RÌ is an element 
of the kernel of IT. Of course, not every f € R™ is the discrete gradient of some x since f may 
contain a curl component. Indeed, a prescribed f€ R™ can be written as a nodal difference 
(f=ITx) if it is cyclically consistent, i.e., if it satisfies Q Tf=0 (i.e., there is no component of the 
flow in the cycle space). Note also that Q is the (cross) differential operator of the graph whose 
expression can be given in terms of a normalized cycle basis; MQ") denotes the subspace of 
R™ with zero flow circulation (curl-free) around cycles. Moreover, I 'f yields nodal accumu- 
lations from flows along edges. As a result, the differential operators, as basic tools of the 
discrete calculus, have been established and properly phrased on the discrete space. This 
mathematical abstraction meaningfully captures the topological structure of the underlying 
discrete setting. Note that the topological characterization of the graph is essentially embodied 
in the algebraic structure of the associated discrete (matrix) operators and their interrelations. 
We also stress the significant distinction between the discrete operators and the commonly 
adopted discretized versions of the continuous differential operations obtained via the method 
of finite differences in numerical analysis; the latter generally lack the desirable topological 
behaviour [42]. 


[-1 1 0 0 1 0 
-1 0 1 0 -1 0 
H=|0 =1 1 0 Q=|1 -1 
0 -1 0 1 0 1 
[0 0 -1 1 0 -1 


Figure 1. An example graph shown along with its edge-node incidence matrix, IT, and cycle (edge- 
Figused inndeneey inarnd hostel erst watsitoedsengqeancidence matrix, IT, and cycle (edge-face incidence) ma 
trix, Q. Note that M=5, Q=4 and R=2. 


As afinal remark $ ne onedd avons oht he Hral inimnod larity, y PEC US property in erent to, the 
when Hee e 1s, inclu in Sack two face that fra hek e in O site 
MAPS ABRSAOES POLAT EAL emely in Rint E SBS PGT In GRis circumstance, t0 efer 
We renallithiahanaty a seAbotidhe ARERR ok every, ROAR SPALA EE: Babe 
stanhalheedaannode ine PACS patriata AARC OR UREA HAP facore {an ipcidleracematrix 
in genatalnis AAA (Har drddow exesnthares gerdacebesidRace maakeratiohUow ga csanherlge 
optimization problem (see also optimization problems (21) and (22) in the following) will be 
integer. Nonetheless, TU property has a further practical significance since the relaxed 
problem, obtained by neglecting the integer constraints, can also be solved using generic 
linear (not integer) programming solvers. 


2.2 Rigorous gradient-based formulation of the MCh-PhU problem 
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is included in exactly two faces that traverse the edge in opposite directions (e.g., a planar 
graph with a minimum cycle basis [30]). In this circumstance, total unimodularity of the edge- 
face incidence matrix stems from the fact that the face-edge incidence matrix is the edge—node 
incidence matrix of the dual graph [30]. Indeed, a TU constraint matrix (and integer constraints) 
guarantees that the solution of the related optimization problem (see also optimization 
problems (21) and (22) in the following) will be integer. Nonetheless, TU property has a further 
practical significance since the relaxed problem, obtained by neglecting the integer constraints, 
can also be solved using generic linear (not integer) programming solvers. 


2.2. Rigorous gradient-based formulation of the MCh-PhU problem 


Once the basic concepts of the discrete calculus and graph theory are presented, we are in the 
position to frame the formulation of the MCh-PhU problem on an appropriate mathematical 
playground. Let us consider a set of Q single-look-complex (SLC) SAR data acquired over a 
certain area of interest. One of them is assumed as the reference (master) image, with respect 
to which the images are properly coregistered. This set is characterized by the corresponding 
acquisition times {f,, ..., fg} and perpendicular baselines {b} ..., bo}. Accordingly, for each 
coregistered SLC pair, a multilook phase interferogram (suitably depurated by the flat-Earth and 
topography contributions, by using a priori information about the topography and acquisition 
geometry) can be produced [37]; however, a common practice (within the multitemporal SB 
InSAR class) is first to identify a suitable small-baseline subset of the relevant multibaseline 
(temporal and spatial-perpendicular baselines) interferometric-pair set [6]. This is done to 
confine the effect of decorrelation phenomena associated with inherent angular and temporal 
electromagnetic backscattering variations [38]. Furthermore, a subset of common pixels of the 
M interferograms are then usually identified via the estimated coherence [37,38], so that only P 
pixels characterized by relatively high coherence values are singled out. In other words, the 
coherence index, providing a quantitative estimation of the decorrelation effects, permits to 
discriminate in favor of the “reliable” pixels. 


The final aim is to reconstruct the absolute (i.e., not restricted in the principal [—7t,7t) interval) 
interferometric phases values from the wrapped (i.e., observed only in the principal [-71,7t) 
interval) interferometric phase pertinent to M multichannel interferograms. 


The problem we are interested in can be naturally framed on a discrete setting. Indeed, one 
possibility is to regard the discrete set of P selected (typically sparsely distributed) coherent 
pixels as a set of nodes, ), in the Euclidean (azimuth, range) plane, and the set of Q SAR 
acquisitions representing a set of nodes, Va, in the Euclidean (temporal-baseline, perpendicular- 
baseline) plane [26]. Accordingly, a formal description of the problem at hand can be given in 
terms of a couple of abstract graphs: G=(Ms;E,) and G=(); Ep) where the corresponding edge 
sets E, and €, have to be properly defined. Accordingly, with Q and M, we denote the 
cardinality of V4 and €,, and with P and N the cardinality of 4 and Ep, respectively. Note also 
that defining E, pertains to the M interferometric pairs selection. Defining a meaningful edge 
set for a collection of nodes is now concerned since different criteria can be adopted to achieve 
it. The dimensionality of the ambient space in which the graph is embedded deserves some 
considerations. In this regard, we recall that a graph is called planar if it can be embedded in 
the plane [30-33]. Note also that a graph is not generally guaranteed to be planar, even if the 
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nodes are embedded in two dimensions. Since planarity is important for the workability of the 
implicated optimization procedure with powerful numerical solvers, a typical option to 
preserve graph planarity is resorting to the Delaunay triangulation in the Euclidean plane for 
establishing the edge set from nodes embedded in two dimensions [26]. Note that such an 
option specifically pertains to the solution strategy [26, 35]; nonetheless, our general formu- 
lation applies as well when different edge structures are adopted. Accordingly, once GY and 
& have been somehow defined, the topological properties inherent to each graph are alge- 
braically captured by the related differential operators, which are summarized in Table 1. 


Symbol Quantity Related meaning 
Q Nodes number of Ga Number of SAR acquisitions 
M Edges Number of Ga Number of interferometric pairs 


Dimension of the cycle space of Ga 


II, M x Q incidence matrix of Ga Discrete-gradient operator 
IF Discrete-divergence operator 
Q, M x R cycle matrix of Ga Discrete-curl operator 

P Nodes number of & Number of selected pixels 
N Edges number of & 

L Dimension of the cycle space of & 

II, N x P incidence matrix of & Discrete-gradient operator 
vey Discrete-divergence operator 
Qg N x L cycle matrix of & Discrete-curl operator 


Table 1. Adopted notation 


First of all, we consider the absolute phase relevant to the multichannel SAR acquisition as a 
node variable pertinent to both the graphs G and & ; by using a matrix representation, this 
information can be conveniently arranged in a Q x P matrix @ as follows: 


Pı 
b =[¢',,....Ẹ"]=] 2 (5) 
Po 


where V p €({1, 2,..., P} pe Re encodes in a vectorial manner the pth node variable relevant 
to the graphs ga ; similarly, Y q €{1, 2, ..., Q} pE R? encodes the gth node variable relevant 
to &. 
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Widely adopted global gradient-based PhU approaches, which have historically been devel- 
oped for the single-channel case, generally consist in three processing steps [1, 29]. First, an 
estimation of the (wrapped) phase gradient is obtained; the estimated phase gradient is then 
suitably corrected (in terms of 27 multiples),and subsequently integrated to attain the 
unwrapped (absolute) phase. 


Within the formulation of the MCh-PhU problem we concern [26], a twofold estimation of the 
discrete gradient field is carried out onto the considered two graphs G, and G&, as discussed 
in the following. The stack of the absolute interferometric phases relevant to the M (vectorized) 
interferograms can be formally represented through a P x M matrix denoted by 


Y 
P =i... p] (6) 
Yp 


wherein the P-dimensional vector p” refers to the absolute phase field pertinent to the mth 
interferogram. Accordingly, W is formally related to the absolute (unwrapped) phase matrix 
@ via the discrete gradient operator I, : 


Py =I Ø (7) 


A 


Note also that 


O TO A" | (8) 


By applying the discrete gradient operator II, to the absolute phase of each interferometric 


pair, we obtain 
I SI ay’ ,.;., Lab” | (9) 
As a result, by using Eq. (7), we get 
I, Y =O DI] (10) 


Second, we consider the (wrapped) interferometric phase that is uniquely defined only in the 
principal value range since it is obtained as the phase of a complex function. Hence, it is 
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convenient to formally introduce the non-injective (modulo-27) wrapping operator 
W:pER-> mod (p +n, 2n)-n €[-7, T). It should be noted that the following trivial identi- 
ties hold: 


W(W(A) + W(B)) = W(A + B) = W(W(A) + B) = W(AtW(B)) a 

W(A)=A+22Z b ve) 
where A and B represent two generic matrices and Z is a suitable integer matrix. Given the 
stack of the unknown absolute (unwrapped) interferometric-phases W, the corresponding 
stack of the wrapped phases W can be conveniently expressed in terms of the wrapped discrete 
gradient of (wrapped) observed phase as follows: 


¥" =[WUIT,W (¢')),..., WHW (p"))] = WUT, W (®)) =W, D) (12) 


where we have exploited Eq. (11b). Note also that the pth column of W", i.e. 
b= WL ‘AW (@?)), can be regarded as an estimate of the absolute-phase discrete gradient on 
the graph ga. It is worth remarking that the observed (multilook) interferometric phase can 


however be corrupted by noise [39-41], which is taken into account by considering an additive 
phase noise term D. Accordingly, by using Eq. (11a), we get 


wT =W(WT,®)+ D) =WU(1,@ + D) (13) 


More specifically, whenever a possible spatial filtering (e.g., conventional multilooking 
followed by a noise-filtering step [42]) is independently applied to each SAR interferometric 
data pair, the resulting term D in Eq. (13) implies that the phase interferograms tp)”, with 
me€{1, 2, ...,M}, are no more fully time consistent (in the sense of [26, 43, 44]). To clarify this 
point, we observe that by using Eq. (13), it turns out that 


W(Q{¥") =W(Q) (I, + D+22Z)) =W(QUI,® + Q{D) =W(aiD) #0 (14) 


where we have used Eq. (11b) with A=11,® + D and noted that Q, Z is an integer matrix and 
QI, =0 (according to Eq. (3)). Eq. (14) reads as “the wrapped discrete curl of the interfero- 


metric phase on &% (i.e., pertinent to the ‘temporal’ domain) is generally different from zero”; 
it formally expresses the (modulo-2r) cyclic inconsistency of the multichannel interferometric 
phase inherent to independently filtered SAR interferograms. Note also that Eq. (14) repre- 
sents, within our framework, the generalization (to a wider class of discrete settings) of the 
“phase triangularity” condition in ref. [44], capturing the underlying structure of the problem 
within a suitable matrix formalism. With reference to the mth interferometric pair, the 
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estimated absolute interferometric-phase gradient on the graph & is then obtained by wrapping 
the discrete gradient of (wrapped) interferometric-phase field: g” =W (II,tp”). Thus, by stacking the 
so-obtained absolute phase gradient estimations, we get the N x M matrix 


G=[g!, g, n eM), where 


G=W(I,¥ ) (15) 
Finally, by substituting Eq. (13) in Eq. (15), we obtain 
G=W(IĻĒ )= w( aw (7,0) +d") (16) 
From Eq. (16), by using Eq. (11b), we get 
G =W(IGÐ I, + ID" +27 TZ) =W(I,®" 1, +D") (17) 


where in the last equality we have noted that IZ is also an integer matrix. It should be 


emphasized that, under the assumption D=o, the equality between Eqs. (10) and (17) holds 
only up to an integer matrix multiplied by 27. 


2.3. The MCh-PhU problem as constrained optimization 


In this Section, the nonlinear inversion MCh-PhU problem is reformulated as a (nonlinear) 
constrained optimization problem. According to the presented general formulation, we 
introduce in the following the MCh-PhU problem as the solution of the following matrix 
equation: 


I, ®' 11, + 11,D' +22 K =G (18) 


where the columns of G represent the interferometric-phase pseudo-gradients estimated from 
the observed phase, and K is an (unknown) N x M integer matrix, whose columns represent 
the corresponding (27t-normalized) corrections to be added to the (wrapped) interferometric- 
phase pseudo-gradients in order to recover the absolute interferometric-phase discrete gradi- 
ents. It is worth noting that the term pseudo-gradient is used here to emphasize that the 
integration of the estimated gradient is path-dependent (non-conservative behavior); the term 
residues [1] is also typically used to connote the inconsistency of the estimated phase gradient. 
As a matter of fact, matrix equation (Eq. (18)) describes an ill-posed problem, in which the data 
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G generally do not constrain sufficiently the problem to get a unique solution. Additional 
suitable constraints and a priori assumption have, thus, to be introduced to solve the problem. 
First, for restoring the cyclic consistency (see Section 2.1) of the estimated pseudo-gradients 
pertinent to the graphs @ and Ga, two corresponding sets of (equality) constraints have to be 
enforced, respectively. More specifically, pre-multiplying both sides of Eq. (18) by Qg and 


taking into account Eq. (3), we obtain 
Q)(G-27K)=0 (19) 


Similarly, by premultiplying both sides of the transposed version of Eq. (18) by Q, and taking 


also into account Eq. (3), we obtain 
O%(G-11,D" - 27 K)"=0 (20) 


Constraints stated by Eq. (19) imply that the columns of G-27 K must lie in the null-space of 
QJ. Since the matrix I], represents a basis to span the null-space of QJ (see Eq. (3)), we may 
then write G-2n K= IX, where X is a new variable. Accordingly, the corrected pseudo- 
gradients stack G-27 K is enforced to be a stack of discrete gradients, which can thus be 
unambiguously integrated. Similarly, Eq. (20) implies [G -MD '-2n K ]'=TIkY. Asa result, the 
two sets of constraints, stated by Eqs. (19) and (20), guarantee that the solution of the problem 
is effective in preserving the cyclic consistency (curl-free) property of the corrected gradients 
pertaining to the graphs G and Ga, respectively. As a matter of fact, the solution of Eq. (18) 
cannot be determined by using the two sets of constraints (Eqs. (19) and (20)) only; thus, the 
inverse problem must be first regularized [45]. The minimum-norm methods search for a 
global solution that minimizes a generalized error-norm associated with an optimality 
criterion, so incorporating prior information about the behavior of the solution [1]. Accord- 
ingly, we resort to a regularization approach using /;-norm minimization in weighted version, 
as a specific case of l -norm general formulation. Formally, the MCh-PhU problem may be then 
formulated as a constrained optimization problem for the field of integer corrections: 


K =argmin |K], (21) 
KeZ 
subject to 


oK" =@"|G-.D"| (27)' 
^ al Les) (22) 
QK = Q;G(27) 
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wherein 


M N 
I = Xem 


m=1n=1 


a 


(23) 


represents the weighted /,-norm [46] of the matrix K, C=[¢jn Ji) xm 


matrix, and Z indicates the field of integer numbers. As far as the existence of an integer 
solution for Eqs. (21) and (22) is concerned, it should be noted that the considerations at the 


denotes a suitable weighting 


end of Section 2.1 apply. Since the first matrix equation in Eq. (22) includes a generally not null 
(unwanted) term Q,'D, its fulfillment deserves further discussion. Although the evaluation of 
W (QD) can be obtained according to Eq. (14), however, a full estimation for Q, D is generally 
not a simple task. Further discussion is provided in Section 3. The solution of the optimization 
problem (Eqs. (21) and (22)) is also referred to as the minimum weighted discontinuity solution 
(in a weighted /,-norm sense) [1, 23]. As a matter of fact, finding the global minimum point of 
the problem stated by Eqs. (21) and (22) for an arbitrary pair of graphs is, in general, a difficult 
task. A suboptimal strategy aimed at solving Eqs. (21) and (22) consists in adopting a two-stage 
approach. This is, in particular, the solution strategy implemented through the extended 
minimum cost flow (EMCF) technique [26], in which the edge structure of each considered graph 
is usually defined via a Delaunay triangulation in the Euclidean plane, to take advantage from 
efficient solvers for minimum cost flow (MCF) problems [47, 49, 50]. We remark that the 
distinctive characteristic of the EMCF approach is the extensive use of the computationally 
efficient MCF method. Moreover, a dual-level parallel model for EMCF has also been proposed 
in refs. [35] and [36]. Moreover, different approaches toward full 3D phase unwrapping have 
recently been proposed in refs. [63] and [64]. 


3. Noise-filtering of multichannel SAR interferograms 


In this Section, we review the basic concepts concerning the filtering of noise that corrupts a 
stack of multitemporal SAR interferograms. First, the noise-filtering operation for single- 
channel multilook interferograms is discussed; subsequently, the general framework of the 
multichannel noise-filtering (MC-NF) approach, which is intimately connected with the 
problem of multichannel phase unwrapping, is described. 


3.1. Decorrelation noise in SAR interferograms 


In order to introduce the problem at hand, let us first consider one single-channel SAR 
interferogram obtained starting from two SAR (synthetic aperture radar) images, namely, i, 
and i, acquired (over the same scene on Earth) at two different times, namely, t; and t, 
respectively. The two SAR images can be represented via two complex-valued signals, say 
i(x, r) and i,(x, r), with x and r denoting the two independent spatial variables (with respect 
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to azimuth and range direction, respectively) in the radar geometry. The two complex signals 
can be expressed as follows [29, 51]: 


An 
; =j T 
i(x,r)=y(x,r)e * +n (xr) 


+6r) 


24 
ad ai 
1,(x,r) = 7,(x,1) e +n, (x,r) 


where ôr is the sensor-to-target slant range difference at time t, with respect to time t,, and 
y(x, r) and y,(x,r) are the corresponding (complex-valued) reflectivity functions of the 
illuminated scene at time ¢, and t, respectively. Furthermore, the two additive (noise) 
contributions n(x, r) and n(x, r) describe random quantities that are included in Eq. (24). As 


a result of these noise terms and of the intrinsic random nature of the two images reflectivity 
functions, when the two SAR images are interfered to form a so-called interferogram, i.e., when 
their phase difference y is extracted, the interferometric phase will be noisy. An important 
parameter influencing the quality of the retrieved interferometric phase is the (complex) cross- 
correlation factor between the two involved SAR images, which is typically defined as 


Efi, -i;] 


Korm 


where p €[0, 1], y E[-n, 7), and the asterisk denotes the conjugate complex value. Notewor- 


= pel” 


(25) 


thy, the cross-correlation factor (Eq. (25)) is a complex-valued term that can be decomposed in 
terms of amplitude p (i.e., p= |x |) and phase y. For interferometric SAR images, x can be 
evaluated by performing spatial averaging (known as multilooking) operations on a statisti- 
cally homogeneous area. Indeed, the symbol E[] in Eq. (25), which is representative of the 
statistical expectation operation [52], can then be replaced by the spatial averaging operation. 
The amplitude factor p, which is known to as coherence, accounts for the similarity between the 
two SAR images, whereas w is the multilook interferometric phase. A value for the coherence 
that approaches zero is representative of an uncorrelated scene, whereas coherence value that 
is close to unity corresponds to a noise-free interferogram. 


There are several causes that are responsible for coherence decrease. As matter of fact, the 
cross-correlation factor in Eq. (25) depends on different noise sources, and it can be conven- 
iently factorized as follows [51]: 


x= Khe ` Kiem ` X spa ` Xaop ` KX mis ` Xvol ` a (26) 
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where 


Xthe 1S the contribution of the thermal noise. 


Xtemp accounts for the effects due to (temporal) changes in the complex-valued reflectivity 
function between the two passages of the radar sensor over the illuminated area. The so- 
called temporal decorrelation is very difficult to be statistically modeled being associated 
to complex modifications of the electromagnetic response of the scene: They can be induced 


by human activities and/or natural causes. 


X spa is the term that takes into account the fact that from one SAR image to another the same 
ground resolution cell is imaged from two slightly different looking angles. The change 
change of the looking angle, in turn, leads to a shift between the range spectra of the two 
SAR images, and accordingly, it causes decorrelation since the range spectra of the two 
interfering SAR images are only partly overlapped. It can be shown that range spectra shift 
depend on the perpendicular baseline of the considered SAR data pair, and there is a limit 
value for the perpendicular baseline (known to as critical baseline) for which the two range 
spectra are completely non overlapped (i.e., the images are definitely uncorrelated one 
another) [29, 51]. 


The term Xa4op takes into account of the so-called Doppler decorrelation effects due to the 


fact that SAR azimuth spectra are centered ona specific frequency (Doppler Centroid). When 
two SAR images with considerably different Doppler Centroid values interfere, a decorre- 
lation noise contribution arises from the imperfect overlapping of the two related azimuth 
spectra. Hence, in the case that the two azimuth spectra are not overlapped at all, we have 


Xdop=9- 
Xmis accounts for possible misregistration between two SAR images. 


Xyo1 accounts for volumetric decorrelation effects [53]. 


The multilook operation, leading to the multilook phase y in Eq. (25), reduces the level of noise 


corrupting interferograms, although this is paid in terms of a reduction of spatial resolution 


of interferograms. Multilook interferometric phase can be described by using a random 


quantity and, accordingly, it can be characterized via the knowledge of its probability density 
function. It has been shown in literature [40, 41, 54] that the probability density function (pdf) of 
an L-multilook interferometric phase (with L being the number of averaging samples in the 


averaging window used for the estimation of the statistical average operation involved in the 


calculation of Eq. (25)) can be given in terms of a Gauss hypergeometric function: 


1 2\L 
r[t+3}a-p W 
avT(L) (1-6?) 


p(w)= [eske (27) 
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where y €[—7, n), B=pcos(p —o), p represents the coherence, ,F | denotes the Gauss hyper- 


geometric function, I'(- ) is the gamma function, and y, represents the expected “true” value of 


the interferometric phase. The peak of the distribution is located at Y = ọọ. 


The pdf in Eq. (27) is sketched in Figure 2a for different values of L and in Figure 2b for different 
values of the p, with ))=0. By observing Figure 2a, it is clear that pdfs become narrower as the 
number of looks L increases (as expected). This finding is extremely important because it 
demonstrates that the interferometric phase may be thought to be corrupted by an additive 
noise random signal, namely, v, that has the same pdf as in Eq. (27) but with a zero-mean 
expected value, i.e., we may assume as valid the following additive model for the interfero- 
metric noise [54]: Y= +v. To further investigate about the statistics of multilook interfero- 
grams, we can observe that the validity of Eq. (27) is only restricted to the [-7,71) interval. 
However, this restriction does not apply when the phase signal is directly derived in the 
complex plane instead of the real plane. In the works of Lopez (2003) [55] and Lopez and Pottier 
(2007) [40], a comprehensive analytical derivation of the noise statistics in the complex plane 
is derived. Nonetheless, the Cramér—Rao bound for the standard deviation of multilook phase 
is given by [59] 


(28) 


that shows that standard deviation depends on the coherence p and multilook factor L. Note 
that the phase standard deviation approaches the limit (Eq. (28)) asymptotically as the number 
of looks increases. 


(a) (b) 


Figure 2. Probability density function of the interferometric phase W [rad]: (a) for different values of 


Figarauntbebabfligokn sit (flack dines dterfedolineiS phhiae line ]1Q-gretitiding, 2reoodnge lime)ervotth 
lgoks(L, (1;—bb)ckdinditfered fing ltedlo£ tine, COrrglatiotineoefficiembpylir(Q, Hbl inep Br ditidrdinhealhes— 


ofthe qelate PAE (Pd—blacks Seti Withee g.5—blue line, 0.7—green line, 0.8—orange line) with L = 
4. 
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3.2. Single-channel noise-filtering approaches for multilook interferograms 


In order to mitigate the effects of decorrelation noise artifacts affecting SAR interferograms, 
several noise-filtering techniques, mostly working on single-channel data, have been proposed 
in literature over the years [42, 54, 56, 57]. As shown in previous Section 3.1, the statistics of 
multilook interferograms can be characterized via a probability density function expressible 
in closed form (Eq. (27)), and the noise standard deviation generally depends upon the 
coherence p and the number of looks L [see also Eq. (28)]. Three different multilook interferograms, 
which are characterized by the same perpendicular baseline (of about 100 m), have been 
obtained by using three SAR sensors working at the different (C, X, and L) bands of the 
microwave region and are depicted in Figure 3. As it is evident from Figure 3, L-band inter- 
ferograms are less affected by noise than the ones pertinent to C and/or X-band. 


Range 


Azimuth Azimuth 


Range 


Azimuth 


Figure 3. Multilook interferogram computed by using different SAR data pairs: (a) July 11, 2011-August 16, 2011, X- 
band Cosmo-SkyMed (CSK); (b) September 15, 2004—October 29, 2004, C-band ASAR/ENVISAT; (c) July 30, 2007-Sep- 
tember 14, 2007, L-band ALOS/PALSAR-1. 
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It should be emphasized that coherence and noise levels can also significantly change from 
one SAR system to another depending on the operational wavelength.Multilook processing 
has been proved to be effective for noise reduction, but this is paid in terms of a decrease of 
the image spatial resolution. Noise filtering constitutes a preliminary step before phase 
unwrapping. Indeed, the multilook operation usually involves an averaging on neighboring 
SAR pixels, hence reducing the spatial resolution of the interferograms. Several algorithms 
have been proposed in literature. The most commonly used noise filter is the boxcar filter 
applied in the complex plane. Another frequently used option is provided by the Golstein’s 
frequency-domain algorithm [42], which is an empirical approach proposed for geophysical 
applications. In this case, a complex interferogram (amplitude and phase) is segmented into 
overlapping rectangular patches and for each patch the relevant power spectrum Z is com- 
puted. 


Range 


Azimuth Azimuth 


Range 


Azimuth Azimuth 


=f Tt 


Figure 4. Multilook interferogram relevant to the SAR data pair July 11, 2011—-August 16, 2011, X-band Cosmo-SkyMed 
(CSK): (a) Original, (b-d) Goldstein filtering, with (b) a =0.25, (c) a@=0.5, and (d) & =1.0. 
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The response of the Golstein’s filter is then computed from the power spectrum as follows: 


Hlé) =| (29) 


where € and 7 denote the relevant spectral variables, respectively. Notice that the filtering 
effect vanishes when a=0; conversely, the filtering effect is more pronounced as the parameter 
a approaches unity. We show in Figure 4 the result of the application of the Goldstein’s filter 
to a multilook interferogram, relevant to the Mt. Etna (Italy) volcano, obtained by using the 
Cosmo-Skymed sensor of the Italian Space Agency (ASI). Specifically, different values of the 
filtering parameter a have been considered in Figure 4. The limited effectiveness of the filtering 
capabilities of Goldstein approach is evident from the result depicted in Figure 4. A modifi- 
cation of the Goldstein filter that relies on an adaptive choice of the filtering factor a (which 
depends on the spatial coherence p) has also (more recently) been proposed by Baran in 2003 
[58]. Other filters, such as the median filter [59] and the two-dimensional Gaussian filter, are 
also used to reduce noise while performing multilooking operations. It is worth noting that 
boxcar and Goldstein filters do not adapt to the direction of the fringes because these filters 
are operated in a square window. In order to overcome such a limitation, Lee et al. 1998 [54] 
then proposed a self-adaptive filter based on local gradient slope estimation that is able to 
improve noise-filtering performance by exploiting directional characteristics of an InSAR 
interferogram. Several adaptations and relevant improvements of the Lee filter have subse- 
quently proposed in literature over the recent years [56, 57], most of them based on the 
exploitation of the intrinsic directional behavior of InSAR interferograms. In fact, compared 
with the fringe phase and gradient, the fringe direction variation is gently, thus making it 
possible to use fringe direction to guide interferogram filtering. 


3.3. The multichannel noise-filtering (MCh-NF) algorithm 


The noise-filtering methods discussed in the previous Section have historically been developed 
to analyse and filter out the noise affecting single interferograms, separately, thus without 
taking into account their mutual temporal relationships. A multichannel noise-filtering 
problem arises when a stack of SAR interferograms need to be jointly exploited. In this case, 
it is profitable to develop/use noise-filtering approaches that not only exploit spatial/frequency 
information but can also take into account temporal relationships among available multichan- 
nel interferograms, in order to identify and filter out the noise affecting the whole interfero- 
metric data stack in the more reliable way as possible. A specific multichannel noise-filtering 
(MCh-NF) method [43], which is based on using a stack of time-redundant multilook inter- 
ferograms, is described in this Section. The MCh-NF method is here described by adopting the 
same rigorous formalism and terminology used for the topological description of multichannel 
phase unwrapping problem presented in Section 2. According to the adopted symbolism, let 
us consider Q SAR images and let M be the number of multilook interferograms characterized 
by small perpendicular and temporal baselines. 
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The resulting interferometric data stack of the M (wrapped) small-baseline multilook inter- 
ferograms can be expressed as p= pi, shay p3] ; thus, the M-dimensional vector p; described 
the (vectorized) multichannel interferometic-phase pertinent to the pth pixel, with 
p E{1, ..., P} and P denoting the number of coherent pixels common to all interferograms. In 
particular, ¥ can be expressed in terms of discrete gradient II,, according to Eq. (13), as: 


PY" =W(IT,® + D) (30) 


wherein ® represents the (unknown) phases associated with the available SAR images, and 
the matrix D describes the additive noise-term that corrupts the stack of interferograms. The 
noise term should be estimated and properly filter out from the generated interferograms. As 
discussed in Section 3.2, the term D arises since both a multilook operation and a noise-filtering 
procedure are typically applied to each single interferogram, separately. Both these operations 
are independently carried out on each single interferogram; hence, they are not necessarily 
time consistent. The fact that the interferograms are not fully time consistent can be formally 


expressed, according to Eq. (14), in terms of discrete curl QJ, in the form: 
W(2;¥")=W(Q{D) #0 (31) 


which represents the topological generalization of the phase-triangularity condition exploited 
by the SqueeSAR technique [44]. Therefore, the multichannel noise-filtering (MCh-NF) approach 
suitably addresses the temporal inconsistencies inherent in the time-redundant multilook 
interferograms, which can be mathematically described in terms of the (modulo-27t) cyclic 
inconsistency of the multichannel interferometric phases [see Eq. (31)]. More specifically, MCh- 
NF is based on the solution of the a nonlinear optimization problem, as detailed in the following. 
First, V p €{1, ..., P}, the Q-dimensional vector (Q is the number of SAR acquisitions) repre- 
senting the (unknown) wrapped phases ®’=W(®’) is estimated as follows: 


2 = (PT wT, db? 
@’ = arg max Ç, oe" pW D )) (32) 


e RÈ 


where j=,/-1 denotes the imaginary unit, P denotes the number of coherent pixels to which 
the noise-filtering procedure is applied, the symbol represents the Hadamard product, and 


C= Lee me, o |’ is an M-dimensional normalized weighting vector representing our confidence 
on the quality of the exploited M (small-baseline) interferometric phases pertinent to the pth 


pixel, with 
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op = (33) 


wherein the generic elements Gr can be related to the spatial coherence as detailed after. 
A 


Subsequently, these estimated vectors D ” are used to reconstruct a new (noise-filtered) stack 
of time-consistent interferograms ¥'=[i], ..., Pp], where T; =W(II,®") and pE{1, ..., P}. 
We emphasize that, according to Eq. (32), the MCh-NF technique is based on searching for the 
(unknown) wrapped-phase vector D”? € RÌ that minimizes the (weighted) circular variance of 
the random (phase) vector representative of the phase difference, pie -W (II,®"), between the 


“original” and the “reconstructed” interferograms. 


The evaluation of the weights for the optimization problem in Eq. (32) is now addressed. Let 
O= [O,, j] be a matrix description for a generic 2-D phase map, whose corresponding vectorized 
representation is provided by the P-dimensional vector ©. Each pixel of the phase map is 
identified by discrete range and azimuth coordinates, denoted by i and j, respectively. 
Accordingly, each pair (i, j) is uniquely associated with a monodimension index p €{1, ..., P} 
identifying an element of the vector O. The spatial coherence relevant to © (i.e., ©) evaluated 
around the pixel (i, j) (associated with the index p) is defined as 


1 Lre La 


PO rager, U] a 


where 2Lp +1 and 2L, +1 are the number of azimuth and range pixels within the used boxcar 
averaging window, which is centred around the generic pixel identified by the discrete range 
and azimuth coordinates, i and j, respectively. In particular, the mth weight Ç," is expressed, 
according to Eq. (34), in terms of spatial coherence relevant to the (vectorized) interferograms 
ip” and evaluated around the pixel associated with the index p, in the functional form 


Cy te) Vme€({l, ..., M}. Therefore, c =G; e, an can be evaluated in terms of the 


spatial coherence directly from the stack of M multilook interferograms W=[h', ..., P~]. As 
experimentally demonstrated in ref. [43], the “reconstructed” interferograms with MCh-NF 
are significantly less affected by noise than the original ones. However, a group of the 
reconstructed interferograms, although limited, can exhibit spatial coherence values smaller 
than the ones relevant to the corresponding original interferograms, thus implying that a 
partial corruption of the spatial coherence occurs during the minimization procedure. In 
particular, it happens in correspondence to interferograms that were originally significantly 
coherent, and this is due to the fact that the strategy in Eq. (32) tends to “inject” coherence from 
very coherent to incoherent interferograms, by exploiting the time redundancy of the small 
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baseline data pairs. Accordingly, in order to also preserve the spatial coherence of the very 
coherent interferograms, a simple nonlinear combination between the original and the 
reconstructed interferograms is carried out, thus further increasing the phase quality of the 
whole set of M reconstructed interferograms. In particular, the two sets of interferograms are 
combined through the following (wrapped) weighted averaging operation: 


wp” = arctan 


Cc” osin(tp”) 40% esin( "| 
G ocos( ip") + g -cos| wb") 


vm e{1,...,M} (35) 


where the symbol represents the Hadamard product, and ¢” and co” are two P-dimensional 
vectors. In particular, C”=[C;", =, Cite LG), =, Gp p» is expressed in terms of the spatial 
coherence relevant to the original multilook interferogram ọ”. Similarly, 
"=, T T y) =, Cp(p”)]' is expressed in terms of the spatial coherence relevant to 
the reconstructed multilook interferogram p". The block diagram of the MCh-NF algorithm is 


depicted in Figure 5. A pertinent pseudo-code for computing the filtered interferometric data 
stack is also presented (Figure 6). 


Note that the exploitation of “conventional” small baseline multilook interferograms is the 
distinctive characteristic of MCh-NF approach with respect to other previous solutions, such 
as the SqueeSAR [44] and Phase Linking [62] methods and other recently proposed multitem- 
poral-filtering techniques [60, 61] based on constraining the analysis to distributed scatter- 
ers [29], which are identified through a pixel-by-pixel selection procedure performed at the 
full resolution complex SAR image spatial grid. Such a selection permits to rely on the 
distributed scattering hypothesis, under which the probability density function (pdf) of the 
complex-valued SAR image may be regarded as being a zero-mean multivariate circular 
normal distribution, and an appropriate maximum likelihood (ML) estimation step of the 
filtered phase values associated to each SAR acquisition is implemented. On the contrary, 
the presented MCh-NF approach focuses on conventional multilook interferograms generat- 
ed without any a priori pixel selection step. Accordingly, in this case, it is not possible to rely 
on the validity of the above-mentioned distributed scattering hypothesis. Therefore, both the 
phase linking [62] and the phase triangulation of the SqueeSAR [44] algorithms require a 
preliminary identification of the statistically homogeneous pixels (SHPs) on the full-resolu- 
tion range-azimuth grid. In particular, in ref. [44], the selection strategy of these pixels is 
based on the application of the Kolmogorov—Smirnov test to carefully select a homogeneous 
statistical population. Clearly, this requires working at the full resolution spatial scale and 
implies the analysis of the amplitude values of the complex SAR image pixels. A more detailed 
comparison among the presented MCh-NF method and the ones provided in ref. [44] can be 
found in ref. [43]. 
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, , Figure 5. MCh-NF block diagram 
Figure 5. MCh-NF block diagram 


4. Expettmental results 
We present in this section some results we obtained by processing a data set consisting of 39 

We peBeimagethisrackt 208, donreresals cobasteid ibe dhis EDTOGA Si senacd Het weeadretsthbetof 39 
SAR2092 gad (Puget 3080 Prenthe au zdliegted Glalnc INESATTE sseaerchidesdea Mycember 
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L’ Atif cane dtn stidan dives banded AEE tate ening sags ENISA AEM 68 Earth- 
quake dskiaeed mdeni ipere Nge Gri RNA Severe 
indutti sha NEE an in Mee SHA RO r OBESA RRUA, we 
retrieved a stack of 338 small baseline differential SAR interferograms with a maximum 
perpendicular baseline of 400 m and a maximum time span of 2000 days [43]. The interfero- 
grams have been computed by performing a complex multilook operation with 4 and 20 
looks in the range and azimuth directions, respectively. For the interferogram generation, we 
used precise satellite orbit information and a three-arcsecond shuttle radar topography mission 


eV Sere were ss ome Fee 


(step 2 Compute matrix j 
for m «1, M do 


{step 2.1: Compute vector } 


for }Opolagical Characterization and Advanced Noise-Filtering Techniques for Phase Unwrapping of... 
Ç ompute using (34) http://dx.doi.org/10.5772/61847 
end for 
(SRTM) unre — using G9) n to remove the topographic phase 
en or . : 
contribu-----—. --------;, ne ec a —--------o--—-— -lave been prefiltered by applying the 
Goldstein’s filter [42]. 
Replace with: 


Algorithm MCh-NF 


Given: Matrix P" = D, pi] =i"... pT 
{step 1: Compute matrix yr =[pr,...apT] } 
for p =1, P do 
{ step 1.1: Compute vectors t, = poé T } 
for m =1, M do 
Compute 6; =¢, (") using (34) 


end for 
for m=1, M do 


Compute a " using (33) 
end for 
Obtain vector ®” using (32) 
Compute vector P =W(IT,®") 


end for 
Given: Matrices ¥", yr and vectors G” e me Coy 
{step 2: Compute matrix ¥=[y’,...,y"] 
for m =1, M do 

{step 2.1: Compute vector 6" =[£",---,E"]" } 

for p =1, P do 

Compute Ê > =, . (ap") using (34) 
end for 


Obtain vector ” using (35) 


end for 


Figure 6. Pseudo code of the MCh-NF algorithm. 


To investigate the performance of the presented noise-filtering approach, we applied the 
nonlinear minimization procedure in Eq. (32) to the stack of the generated (original) multilook 
small baseline interferograms. As a result, we retrieved a new set of noise-filtered interfero- 
grams that are characterized by generally improved coherence values. This is clearly visible 
in Figure 7a-f, where we compare three unfiltered (left side) interferograms with the corre- 
sponding (right side) noise-filtered interferograms. It is evident how the fringes due to the 
earthquake are well recovered. Such interferometric data stacks can then be used for the 
generation of surface deformation time series using the small baseline subset (SBAS) [6] 
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Figure 7. Comparison between the original (left column) and noise-filtered (right column) multilook interferograms 
relevant to the area of the Abruzzi region (Italy). (a-c) October 30, 2005, to November 8, 2009, August 21, 2005, to June 
6, 2010, and August 1, 2004, to April 12, 2009, interferograms, characterized by perpendicular baseline values of 192, 
145, and 395 m, respectively. (d-f) Noise-filtered multilook interferograms corresponding to the ones in panels a-c, re- 
spectively. 


processing chain, whose parallel version (P-SBAS) has been proposed in refs. [13, 14, 35, 36]. 
This is matter for the analysis presented in the next subsection. 
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small baseline Advanced EMCF-SBAS processing chain 


Figure 8. Block diagram of the advanced EMCF-SBAS processing chain. 


4.1. The use of MCh-NF algorithm within MCh-PhU framework 


We present in this subsection how the MCh-NF algorithm can be efficiently used within the 
SBAS processing chain, where phase unwrapping procedures are implemented through the 
MCh-PhU technique known as extended minimum cost flow (EMCF), also discussed in the first 
part of the chapter. Figure 8 shows the diagram block of the advanced EMCF-based SBAS 
processing chain [26, 35, 43], which integrates the conventional SBAS codes, exploiting the 
EMCF MCh-PhU procedure to perform phase unwrapping operations, with the presented 
MCh-NF noise-filtering technique. In addition, an effective procedure for the selection of time- 
redundant interferograms is also included; interested readers can find additional details in ref. 
[43]. To provide an example of the potential of the advanced processing chain incorporating 
both the discussed MCh-PhU and MC-NF techniques, we here focus on the area of Yellow- 
stone caldera, representing one of the largest and most active volcanic systems in the world. 
We analyze the temporal evolution of the surface deformation occurring in this area by 
applying the implemented EMCF-SBAS processing chain to a set of 22 ENVISAT images (Track 
41, Frame 2709), acquired from May 2005 to September 2010, from which we have retrieved a 
corresponding set of 122 small baseline interferograms [43]. As in the previous case study 
(relevant to Abruzzi area), the prescribed maximum values of 400 m and 2000 days for 
perpendicular and temporal baseline, respectively, have also been considered. The retrieved 
mean deformation velocity map depicted in Figure 9, which has been obtained by applying 
the processing chain (MCh-NF + EMCF-SBAS) of Figure 8, allows us to recognize the complex 
deformation scenario affecting the Yellowstone Caldera region and its surroundings, where 
uplift effects and broad subsidence patterns characterize the detected deformation field. In 
addition, the deformation time series relevant to four pixels, whose locations are identified by 
the black stars labeled as P1, P2, P3, and P4, are also illustrated in Figure 9. 


5. Conclusion 


Within the context of SAR interferometry, two main issues related to the multichannel phase 
unwrapping and noise filtering for interferometric data stacks processing have been ad- 
dressed. First, a rigorous gradient-based formulation for the multichannel phase unwrapping 
(MCh-PhU) problem has been systematically established, thus providing a topological 
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Figure 9. (a) Mean deformation velocity map (in color) of Yellowstone Caldera, computed in coherent pixels only and 
superimposed on the SAR amplitude image (gray-scale representation) of the zone, retrieved by applying the ad- 
vanced EMCF-SBAS processing chain. The black square marks the location of the reference SAR pixel. (b-e) Deforma- 
tion time series relevant to the four pixels identified via black stars in Fig. 8(a). 


characterization of the problem within the purview of the theoretical foundation of the discrete 
calculus. Then the innovative MCh-NF procedure for the noise filtering of time-redundant 
multichannel multilook interferograms has been properly presented within the considered 
topological framework, by adopting a consistent formalism. Finally, some experimental results 
obtained with real data have been shown, thus demonstrating the effectiveness of our 
approaches and their relevance for geospatial phenomena analysis and understanding. 
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